/*!  \file  commonSVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/03/2013 12:57:34 PM
 *   \copyright GNU Public License.
 */
#include <iostream>
#include <fstream>
#include <assert.h>
#include <float.h>
#include <stdlib.h>
#include "Matrix.h"
using namespace std;
bool checkSVDDimentions(const Matrix & A, const Matrix & U,
				 						 const Matrix & S, const Matrix & VT){
	bool dimentionMatch =   A.rows() ==  U.rows() &&   U.rows() == U.cols() && 
	U.cols() == S.rows() && S.cols() == VT.rows() &&  VT.cols() == A.cols();
	if (!dimentionMatch) 
		cerr << "[void checkSVDDimentions] Error dimentions don't match" << endl;
	return dimentionMatch;
}
bool checkSVDMemory(const Matrix & A, const Matrix & U, 
									const Matrix & S, const Matrix & VT){
	bool haveMem = A.getData() != NULL &&  S.getData() != NULL && 
						U.getData() != NULL && VT.getData() != NULL;
	if (!haveMem) 
		cerr << "[checkSVDMemory] Error Matrices are not allocated" << endl;
	return haveMem;
}
void getResidual(double & resID, Matrix & A, Matrix & U, Matrix & S, Matrix & V) {
   assert(V.isOrthogonal());
   assert(U.isOrthogonal());
	double eps = DBL_EPSILON;
	int n = A.minDim(); 
	double r1 = (A * V.T() - U * S).nrm2();
	double r2 = (A.T() * U - ( S * V).T()).nrm2();
   resID = fmax(r1,r2)/(n * A.nrm2() * eps) ;
	if (resID != resID ) { // is NaN
		cout << A << endl;
		cout << S << endl;
		cout << U << endl;
		cout << V << endl;
		cout << U * S * V.T() << endl;
		exit( -1 );
	}
}
void getOrth(double & resID, Matrix & U) {
	resID = U.OrthNorm() / (U.minDim() * DBL_EPSILON);
}

void getNorms(double * norm, Matrix & A, Matrix & U, Matrix & S, Matrix & V){
	norm[0] =  (        A * V.T() - U * S    ).nrm2();
	norm[1] =  (        A         - U * S * V).nrm2();
	norm[2] =  (U.T() * A         -     S * V).nrm2();
	norm[3] =  (U.T() * A * V.T() -     S    ).nrm2();
	norm[4] =  U.OrthNorm();									
	norm[5] =  V.OrthNorm();									
}
